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Abstract. 

The effects of locally random magnetic fields are considered in a nonequilibrium 
Ising model defined on a square lattice with nearest-neighbors interactions. In order 
to generate the random magnetic fields, we have considered random variables {h} that 
change randomly with time according to a double-gaussian probability distribution, 
which consists of two single gaussian distributions, centered at +ho and —ho-, with the 
same width a. This distribution is very general, and can recover in appropriate limits 
the bimodal distribution {a 0) and the single gaussian one {ho = 0). We performed 
Monte Carlo simulations in lattices with linear sizes in the range L — 32 — 512. The 
system exhibits ferromagnetic and paramagnetic steady states. Our results suggest the 
occurence of first-order phase transitions between the above-mentioned phases at low 
temperatures and large random- field intensities ho, for some small values of the width 
a. By means of finite size scaling, we estimate the critical exponents in the low-field 
region, where we have continuous phase transitions. In addition, we show a sketch of 
the phase diagram of the model for some values of a. 



PACS numbers: 05.10.Ln, 05.50.+q, 64.60.De, 75.10Hk, 75.40.Mg 
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1. Introduction 

The Random Field Ising Model (RFIM) is one of the most studied systems in magnetism 
(for reviews, see [H |2] and more recently [3]). Its simple theoretical definition and 
the interesting physical properties that emerge from its study represent two great 
motivations for the investigation of this model. In addition, a considerable experimental 
interest arised after the identification of the RFIM with some diluted antiferromagnets 
in the presence of a uniform magnetic field |H [5], like FexZni_xF2 and FexMgi_xCl2 



Following the success of the Ising model to capture the essential physics of complex 
systems, several lattice models have been proposed [8]. These models involve a periodic 
lattice whose sities are occupied by spin- 1/2 variables with two states, = ±1. Every 
configuration, s = {sj}, has a potential energy given by the Ising Hamiltonian, 



i.e., any pair of nearest-neighbors spins contributes with a random exchange interaction 
Jij and each spin interacts also with a magnetic field of local intensity h^. Systems based 
on Eq. ([T]) are called models of quenched disorder (MQD) [8J. They are characterized 
by a frozen- in spatial distribution of disorder, i.e., Jij and/or hi vary at random with 
i but remains fixed with time. These models have been very well explored in the 
literature [2], and some experimental results can theoretically be reproduced even at 
mean-field level [9]. However, MQD are equilibrium models, and they neglect diffusion 
of magnetically active ions that occurs in some systems, like spin glasses [2] . Diffusion 
constantly modifies the distance between each specific pair of spin ions in substances 
and, consequently, one should probably allow for variations both in space and time 
of Jij and/or hi in a model. These effects do not seem to be correctly described by 
annealed systems, where the change with time of the spatial distribution of Jij and/or 
hi is constrained by the need to reach equilibrium with the other degrees of freedom 
[T0| [m IT2I \T3\ IT4] . Therefore, impurities tend to be strongly correlated, which is not 
observed in most substances. Instead, one may conceive a situation in which both s 
and the spatial distribution of impurities Jij, or the spatial distribution of the local 
fields hi, vary with time. That is, one may assume that spins and impurities (or fields) 
behave more independently of each other than in the annealed model so that a conflict 
occurs, and a steady nonequilibrium condition prevails asymptotically. This is consistent 
with the reported observation of nonequilibrium effects, for exemple, the influence of 
the details of the dynamical process on the steady state in some materials [8]. These 
models can be called nonequilibrium random-field Ising models (NRFIM) 

The study of nonequilibrium models defined by Eq. ([1]) reveals many interesting 
features [8], with a rich variety of phase transitions and critical phenomena. Exact 
results obtained for one-dimensional lattices [151 [HI [H] Monte Carlo (MC) 

simulations [HI [19] reveal that the critical behavior is non-universal, but it generally 
depends on apparently irrelevant details of dynamics, like diffusion of impurities, i.e., 



mum- 




(1) 



<i,j> 



First- order phase transition in a 2D random-field Ising model with conflicting dynamics?) 



the properties of the distribution of the random variables, and on the transition rates 
chosen. The probabihty distributions analyzed were discrete, for random exchange 
interactions and/or random magnetic fields [8]. As experimental realizations of these 
systems, we may imagine : i) a magnetic material under the action of a random (or 
very rapidly fluctuating) magnetic field, i.e., a field that varies according to a given 
probability distribution with a period shorter than the mean time between successive 
transitions that modify the spin configuration, or ii) a disordered system with fast and 
random diffusion of impurities in which the latter consist of both exchanges as in spin 
glasses and local fields as in random field systems [T5] . 

Due to these motivations, we have studied the NRFIM on a square lattice with 
nearest-neighbors interactions and in the presence of random magnetic fields that follow 
a double-gaussian probability distribution. We performed Monte Carlo simulations and 
our results suggest that first-order phase transitions occur in the model for some values 
of 0" > at low temperatures and high random-field intensities. In the low-field region, 
we have performed a finite size scaling (FSS) to estimate the critical exponents, in order 
to test whether the system follows continuous FSS laws. 



2. Model and Monte Carlo Simulation 



A simple realization of the ideas discussed in the introduction is the Kinetic Ising Model 
where the spin system is in contact with a heat bath at temperature T that induces 
stochastic changes of s according the master equation Ui 

dPt{s) 



dt 



^[c(s^^)Pi(sO-c(s;z)P,(s)]. 



(2) 



Here, -Pt(s) is the probability of any configuration s at time t, and c(s; i) is the 
probability per unit time for a transition from s to s*; the latter is obtained from s 
by flipping spin Sj. The detailed balance condition, i.e., c(s; i) = c(s*;i) exp[— /3A7i], 
with A7i = 7Y(s*) — 7i(s), where (3 = {kT)~^ is the inverse temperature {k is the 
Boltzmann constant), is sufficient to guarantee that the stationary solution of Eq. ([2]) 
is the Gibbs state corresponding to energy given by Eq. ([1]) and temperature T [H [18] . 
This is satisfied by the Metropohs algorithm [20j, c(s; i) = min{l, exp(— /5A7i)}, for 
instance [T8] . 

We have considered the Hamiltonian given by Eq. ([T]) on a square lattice of linear 

J > 0. In other words, 



dimension L and constant exchange interaction, i.e., Jij 



(3) 



<i,j> i 

and the random fields {hi} follow a continuous probability distribution, namely, the 
double-gaussian one. 



2V2 



exp 



jh^ - hof 
2a2 



+ exp 



{hi + hp) 
2ct2 



(4) 
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where we have two gaussian distributions centered at ±/io with the same width a. The 
time evolution of any state s is generated by Eq. ([2]) with c(s; x) describing a competing 
process [18] . For the numerical implementation, the algorithm is as follows: at each time 
step, a new configuration of random fields {hi} is generated according to P(/ij), Eq. (jl]); 
then, every lattice site is visited, and a spin flip occurs according to Metropolis' rule. 
In other words, the random variables hi vary with time, i.e., the system is described 
at each time by Eq. ([3]), with hi distributed according to P{hi) given by Eq. (jl]). 
The instantaneous Hamiltonian, Eq. (13]), corresponds to the one that characterizes the 
quenched RFIM. However, the local fields also change at random with time at each 
site according to P{hi). Consequently, a sort of dynamic conflict occurs, which differs 
fundamentally from the equilibrium quenched and annealed random-field models, and 
we have a nonequilibrium system [8]. We may interpret this as a magnet under a field 
that is changing (all over system) at random with time according to P{hi) [15] [T6] . 
Thus, we have two different characteristic time scales: one for the fluctuations of the 
spins and another one for the fluctuations of the random field, and we consider that 
these two fluctuations are independent (a formal discussion about this is found in ref. 
[H], chapter 7). We will show that this rapid fluctuation of the random variables hi 
during the time evolution of the system lead to a flrst-order transition between the 
ordered and the disordered phases for high values of the random-fleld intensity ho, in 
opposition of the continuous transition that was found for small ho. 

Although we do not have evidences that the double-gaussian distribution for the 
random flelds can be realized experimentally, we can argue that this distribution is 
suitable for an appropriate theoretical description of random-fleld models [9]. In the 
identiflcation of the RFIM with diluted antiferromagnets in the presence of a uniform 
magnetic fleld, the local random flelds are expressed in terms of quantities that vary in 
both signal and magnitude [U [5] . This characteristic rules out the bimodal probability 
distribution from a such class of physical systems. Although the RFIM deflned in terms 
of a single gaussian probability distribution is physically acceptable, it usually leads only 
to continuous phase transitions, either within mean-fleld [21 1 [22 ] 123]. or standard short- 
range- interaction approaches [21] [25]. This discussion is based on quenched random-fleld 
models, but an extension to the case of nonequilibrium diluted antiferromagnets was 
discussed by some authors [H][lH]- As pointed in a previous paper [9], the double-gaussian 
distribution can recover in appropriate limits the bimodal distribution {a — > 0) and the 
single gaussian one {ho = or a > ho). In Fig. [T] we show two histograms of the random 
variables {hi} generated numerically accordingly the probability distribution of Eq. (jl]). 

In the following we use for simplicity J = 1. We have studied systems of 
L = 32, 64, 128, 256 and 512 with periodic boundary conditions and a random initial 
conflguration of the spins. We have analyzed the following values of the parameters: 
a = 0.0,0.1,0.2,0.3,0.4 and 0.5, all of them in the range 0.0 < ho < 3.0 and 
0.01 < T < 4.0. The results for all values of the magnetic fleld parameters a and ho show 
that flnite-size effects are less-pronounced for L > 128. We can test the equilibration 
of the system by monitoring the magnetization as a function of the MC time. We have 
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Figure 2. Evolution of the magnetization as a function of the Monte Carlo time, for 
L — 128 and a = 0.2. Examples for low (left side) and high random field intensity 
ho (right side) are shown. As discussed in the text, the equilibrium states are easily 
anchieved for low ho but arc harder to obtain as ho increases. 



found that the equihbrium is easily achieved for low ho but is harder to obtain as the 
intensity of the random field increases [see Fig. [2]. Thus, we have used 10^ MC steps 
for equilibration and 10^ MC steps for averaging. In addition, we have to consider the 
autocorrelation function, defined as [261 



< m{t + to)m{to) > - < m >^ 
< > - < m >^ 

where m{t) is the value of the magnetization of the system at time t. The averages 
in Eq. ([5]) are over all possible time origins to for a steady state. For sufficiently 
large t, m(t) and m(0) will become uncorrelated, and hence < m{t + to)m{to) >— >< 
m{t + to) >< m(to) >=< m >^, i.e., C(t) — > 0. In the above-described dynamics, the 
autocorrelation function C{t) decays in the traditional exponential form, C{t) ~ e"*/"^, 
where r is the autocorrelation time, whose magnitude depends on the choice of the 



First- order phase transition in a 2D random-field Ising model with conflicting dynamics^ 



1" 



0.8 
0.6 

a 

0.4 



0.2 



"0 1 2 3 4 

T 

Figure 3. Magnetization per spin as a function of temperature for L — 128, a = 0.0 
and some typical values of the magnetic field intensity hg. For high values of ho the 
system is in the paramagnetic phase, but for all other values of ho we can observe 
a continuous transition between the ordered and the disordered phases. The symbols 
(circle, square, diamond, star and triangles up and down) are computed points, whereas 
the lines are just guides to the eye. The error bars are not shown because they are 
smaller than data points. 




physical parameters of the system (T, /iq, a). Thus, to obtain configurations that are 
statistically independent, we have used 30 MC steps between two measurements of the 
quantities of interest (remember that in the standard 2D Ising model, 10 MC steps are 
sufficient). 

3. Results 

In the following we will show results for the order parameter m, the magnetization per 
spin, and for x i the magnetic susceptibility, which can be obtained from the simulations 
by the fluctuation-dissipation relation, 

< > — < m >^ 
X = ^ , (6) 

where k is the Boltzmann constant and < > stands for MC average p^. For all following 
results we use k = 1. 

Our flrst simulations were done for a = 0.0, the limit case of the bimodal 
distribution. In Fig. [3] we show the results of the magnetization per spin, m, versus 
the temperature, T, for L = 128 and some typical values of ho- These results show that 
the system is in the disordered paramagnetic phase for high values of ho (> 2.0), but 
for all values ho < 2.0 there are continuous phase transitions between the ferromagnetic 
(m 7^ 0) and the paramagnetic (m = 0) phases. 

We can analyze now the behavior of the magnetization for a > 0.0. In Fig. H] we 
show the magnetization per spin as a function of the temperature for L = 128, a = 0.2 
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Figure 4. Magnetization per spin as a function of temperature for L = 128, a = 0.2 
and some typical values of the magnetic field intensity ho- We can observe that m 
decreases continuously with T for small values of /iq, but we can identify jumps on the 
magnetization for high magnetic-field intensities. The error bars are smaller than data 
points. 
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Figure 5. Magnetization per spin as a function of temperature for a = 0.2, ho = 2.0 
and different lattice sizes. We can observe jumps on the magnetization for all lattice 
sizes. The error bars are smaller than data points. 



and some typical values of h^. In this case, for small ho we can observe the same behavior 
as in the case of a = 0.0. The difference appears for higher values of the magnetic-field 
intensity: we can observe jumps in the magnetization for values near ho = 2.0. So, the 
next step is to verify the behavior of the system for a = 0.2, ho = 2.0 and different linear 
lattice sizes L. In Fig. [5] we show the results of magnetization versus the temperature 
for L = 32, 64, 128, 256 and 512. As discussed before, for the smaller value of L we 
have pronounced finite-size effects. Thus, for the calculations of the exponents, we will 



First- order phase transition in a 2D random-field Ising model with conflicting dynamicsS 




t ' log L 



Figure 6. In the upper figures we sliow tlie susceptibility % as a function of the 
temperature for a — 0.2, ho — 2.0 and L = 128 (left side) and the susceptibility peaks 
Xmax versus L in the log- log scale for some lattice sizes (right side). In the lower 
figures we also show the susceptibility x as a function of the temperature for a — 0.2 
and L ^ 128, but for ho — 0.5 (left side) and the susceptibility peaks Xmax versus L in 
the log-log scale for some lattice sizes (right side). The straight lines have slope equal 
to 2.087 (for ho = 2.0) and 1.7157 (for ho — 0.5), respectively. As explained in the 
text, the peaks grow fast and we only show the result for x for one lattice size. We 
have discarded the smaller lattice size L ~ 32 for the first-order transition calculation 
(upper figure, right side), as explained in the text. The error bars are smaller than 
data points. 

discard the size L = 32. 

The jumps observed in the magnetization in Fig. [5] suggest a first-order 
ferromagnetic-paramagnetic phase transition at low temperatures. In order to verify 
this suggestion, we shall analyze the susceptibility X; Eq. ([6]). We have verified that 
the susceptibility peak grows very fast with the linear lattice size L, which rules out 
the possibility of the plot of x versus T for all lattice sizes studied. Due to this, we 
have plotted in Fig. [6] (upper figure, left side) the susceptibility as a function of the 
temperature only for the size L = 128. Also in Fig. [6] (upper figure, right side) we show 
the log-log plot of data points of Xmax versus L for L = 64, 128, 256 and 512 (we have 
excluded the smaller size L = 32 due to finite-size effects, as observed earlier). If a 
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Figure 7. Binder cumulant C/l, Eq. ([9]), for a = 0.2, ho = 0.5 and some lattice sizes 
as a function of temperature (left side) and the best collapse of data (right side), based 
on Eq. (fTO|) . The critical values are Tc = 2.1258 and v — 1.0277, as discussed in the 
text. The error bars are smaller than data points. 



first-order transition occurs in the system, the peaks Xmax of the susceptibility scale as 

Xmax ~ L\ (7) 

where d is the dimension of the lattice [27] (in our case, d = 2). Fitting data, we have 
observed that Xmax scales as Xmax ~ L°-i where 

a = 2.087 ±0.036, (8) 

which is compatible with the form given in Eq. ([71). Hence, for a = 0.2 and ho = 2.0, 
our simulations suggest that a first-order phase transition occurs in the system for low 
temperatures. We have verified that a similar behavior (jumps in the magnetization, 
with susceptibility diverging with L^) also occurs for other values of the width of the 
distribution a, like a = 0.1 and 0.3 (we do not show these results because they are 
similar to the magnetization curves of Fig. [5|) . Numerically, it is difficult to determine 
the threshold value acrit for which we only have continuous phase transitions, as it was 
done analytically in [9], for the equilibrium (quenched) case. However, for all of the 
simulations with a > 0.3 we did not detect first-order transition features. 

We have performed a preliminary calculation of the critical exponents in the low- 
field region in order to test whether the system follows continuous FSS laws. Thus, 
we have choosen the values a = 0.2 and ho = 0.5, for which our simulations suggest a 
continuous phase transition (see Fig. HI), and we have calculated the magnetization and 
the susceptibility for various lattice sizes. The critical temperature of the infinite lattice 
Tc = Tc(oo) was obtained by extrapolating the Tc{L) values given by the susceptibility 
peaks positions, and for this case we have = 2.1258 ± 0.0015. The exponent related 
to the divergence of the correlation lenght u was calculated by means of the Binder 
cumulant [28], defined as 
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Figure 8. Magnetization per spin as a function of temperature (left) and the scaling 
plot of the magnetization (right). The best collapse was obtained for Tc = 2.1258, 
/3 — 0.0688 and ^ — 1.0277, as discussed in the text. The error bars are smaller than 
data points. 



which has the FSS form 

Ul = Ul{{T-T,)L^I^), (10) 

where Ul is a scahng function that is independent of L. In Fig. [7] we show the Binder 
cumulant for some lattice sizes (left side) and the best colapse of data (right side), based 
on Eq. ffTOj) . obtained with the above-mentioned value of T^. and v = 1.0277 ± 0.0102. 
The exponent f3, which characterizes the behavior of the magnetization near the critical 
point Tc, was calculated to give us the best collapse of the magnetization curves. In 
Fig. M we show the magnetization per spin versus the temperature (left side) and the 
best collapse of the curves (right side), based on the standard FSS forms 

T,(L) = T, - a 

m{T,L) = L-f^/TiiT -T,) L^/"), (11) 

where a is a constant and F{x) is a scaling function that has the limiting behaviors 
F{x << 1) ~ and F{x » 1) = c, where c is a constant. From the above FSS forms, 
Eqs. ( fTTI) . we have found that P = 0.0688 ± 0.0032. We have verified again that the 
susceptibilty peak grows very fast, and we show in Fig. [6] (lower figure, left side) the 
susceptibility as a function of the temperature only for L = 128, in order to compare 
the shape of this curve with the one obtained for a first-order transition (Fig. [6l upper 
figure, left side). Also in Fig. [6] (lower figure, right side), we show the peaks of the 
susceptibility versus the linear lattice size L in the log-log scale for L = 32, 64, 128, 256 
and 512. For the continuous transition, the finite size effects are not so pronounced as 
in the first-order case, and we have considered even the smaller lattice size, L = 32. 
These peaks scale with the system size in the traditional power-law form for continuous 
transitions, i.e., 



Xmax ~ L 
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Figure 9. Sketch of the phase diagram of the model, separating the paramagnetic (P) 
and ferromagnetic (F) phases, for a = 0.0, 0.2 and 0.5. The symbols (circle, square 
and diamond) are computed points, whereas the lines are just guides to the eye. Full 
(dotted) lines represent continuous (first-order) phase transitions. As discussed in the 
text, we have first-order phase transitions for low temperatures and high random- field 
intensities only for some small values of cr. In the limit of the bimodal probability 
distribution (cr 0) and for large cr, we only have continuous phase transitions. The 
error bars are smaller than data points. 



where 7/z/ = 1.7157 ± 0.0499. For the value of u obtained above, we have 7 = 
1.7632 ± 0.0005. In other words, the standard FSS forms for continuous transitions, 
Eqs. ffTOj) and ffTTl) . are the correct FSS equations for that nonequihbrium model. It 
is the purpose of this exponents' calculation; how the critical exponents depend on the 
choice of the parameters of the disorder distribution (cr and ho) is a question beyond 
the target of this work. To summarize the results discussed along this paper, we show 
in Fig. [9] a sketch of the phase diagram of the model, in the plane temperature T versus 
random field intensity ho for some values of the parameter a. 

4. Conclusions 

In this work, we have studied the nonequilibirum random-field Ising model on a square 
lattice with nearest-neighbors interactions by means of Monte Carlo simulations. The 
lattice sizes analyzed were L = 32, 64, 128, 256 and 512. Due to computational cost, we 
have used only 10^ MC steps per spin for averaging. In order to generate the random 
fields, we have used the recently proposal double-gaussian probability distribution, which 
consists of two single gaussian distributions, centered at -\-ho and —ho, with the same 
width cr. As pointed in a previous work |9], this distribution is expected to be more 
appropriate for the theoretical description of real systems than other simpler well-known 
cases, i.e., the bimodal and single gaussian distributions. 
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The time evolution of the system is stochastic because of a competing spin-flip 
kinetics, which, in addition to the usual heat bath, involves a random external magnetic 
field. The competition induces a kind of dynamical disorder that might be present in 
real disordered systems such as the class of random- field materials [16j. This system 
differs from the standard equilibrium ones: while the local field is randomly assigned in 
space according to a distribution P{hi), which remains frozen in for the quenched case, 
and P{hi) contains essential correlations in the annealed system, where the impurity 
distribution is in equilibrium with the spin system, our case is similar to the quenched 
system at each time during the stationary regime, but hi keeps randomly changing with 
time, also according to P{hi), at each site i. 

Our results suggest that in the limit of the bimodal distribution [a = 0.0), the 
transition from the ordered to the disordered phase is continuous. However, first-order 
phase transitions may occur in the system for low temperatures and high random-field 
intensities ho, for some values of the parameter a. The order of the transition was 
determined by the scaling of the susceptibility peaks, which must grow with the total 
number of spins in the case of first-order transitions j27]. The threshold value of a, for 
which the phase transitions are always continuous, is difficult to determine numerically. 
Nonetheless, the simulations suggest that for a > 0.3 first-order transitions do not occur 
anymore. 

We have performed a preliminary calculation of the critical exponents in the low- 
field region in order to test whether the system follows continuous FSS laws. So, we 
have choosen the values a = 0.2 and ho = 0.5, for which our simulations suggest a 
continuous phase transition, and we have calculated the critical exponents. We have 
found that the system obeys standard continuous FSS laws, with /3 = 0.0688 ± 0.0032, 
u = 1.0277 ± 0.0102 and 7 = 1.7632 ± 0.0005. Note that the susceptibility peaks grow 
with the system size L with an exponent 7/z/ < 2, as expected for continuous phase 
transitions in two-dimensional lattices. 

The mean field analysis of the Random Field Ising Spin Glass (RFISG) that 
follows the double-gaussian probability distribution was done recently for the case 
of quenched variables That study showed a rich behavior of the system, with 

continuous and first-order phase transitions as well as a change in the concavity of 
the Almeida-Thouless (AT) line [30], which was experimentally verified in the diluted 
antiferromagnets FexZni_xF2 [31]. However, recent numerical results suggest that there 
is no AT line for short-range models with bond disorder and/or under random magnetic 
fields [32l |33| Thus, extensions of this work concerning on MC simulations of the 
RFISG in the presence of random fields generated according to the double-gaussian 
distribution, for the equilibrium (quenched and annealed variables) and nonequilibrium 
cases, would be of great interest to elucidate the existence of the AT line in equilibrium 
systems with short-range interactions and to obtain more results for nonequilibrium spin 
glasses. 
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